Capillary Waves at Liquid/ Vapor Interfaces: A Molecular Dynamics Simulation 
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Evidence for capillary waves at a liquid/vapor interface are presented from extensive molecular 
dynamics simulations of a system containing up to 1.24 million Lennard- Jones particles. Careful 
measurements show that the total interfacial width depends logarithmically on Lii , the length of the 
simulation cell parallel to the interface, as predicted theoretically. The strength of the divergence of 
the interfacial width on Lu depends inversely on the surface tension 7. This allows us to measure 7 
two ways since 7 can also be obtained from the difference in the pressure parallel and perpendicular 
to the interface. These two independent measures of 7 agree provided that the interfacial order 
parameter profile is fit to an error function and not a hyperbolic tangent, as often assumed. We 
explore why these two common fitting functions give different results for 7. 

PACS number(s): 68.35.Ja 68.35.Md 64.70. Fx 68.35.Ct 68.10.-m 



An interface is the physical boundary between two dis- 
tinct thermodynamic phases, i.e. a region characterized 
by a local gradient of the order-parameter which mean 
value changes from one phase to the other. Examples in- 
clude domain boundaries in ferromagnetic materials, the 
interface between two immiscible liquids, or between a 
liquid and its own vapor below the critical temperature 
T c . This last case has been well studied, both theoret- 
ically and experimentally. For simple fluids interacting 
via van der Waals forces, the mean local density changes 
monotonically [jjj across the interface from its bulk liquid 
value to that of the vapor. In other systems, such as alkali 
metals for example (|,|], the profile across the interface is 
often more complex, with oscillations in the local density 
superimposed on the decaying density profile. 

For simple fluids, thermodynamic considerations alone 
would predict that the interfacial width w, depends only 
on temperature and on the interaction energies within 
each phase and across the interface. However, the pres- 
ence of the interface breaks the translational invariance 
of the system, inducing Goldstone fluctuations or "cap- 
illary waves" at an interface [f|j|]. For two-dimensional 
interfaces, these non-critical fluctuations give rise to a 
logarithmic increase in the interfacial width w with in- 
creasing Lu , the length of the interface. Evidence for cap- 
illary waves has been found experimentally from X-ray 
scattering p- p[ o n liquid/ vapor interfaces and neutron 
reflectivity on polymer/polymer interfaces. More- 

over, nuclear reaction analysis (NRA) depth profiling |l2] ] 
has been used to directly investigate the film thickness 
dependence on the interface width between two polymer 
films and is in qualitative agreement with capillary-wave 
predictions. Capillary waves have also also been observed 
in computer simulations for polymer/polymer interfaces 
[H-|l6|. Most previous simulations p|]T^| of the liq- 
uid/vapor interface in three dimensions did not inves- 
tigate the dependence of w on the size of the interface. 
One recent simulation study fl9| of a thin polymer-film 



system gave some evidence for capillary waves, but the 
longitudinal size of the interface was very small. 

The purpose of this paper is to present computer sim- 
ulation results of interfaces in a liquid/vapor system. To 
our knowledge, these simulations are the most exten- 
sive studies of the interface fluctuations due to capillary 
waves. In particular, we obtain the surface tension 7 two 
different ways: from the dependence of w on L\\ (7™), and 
from the difference in pressure parallel pii and perpendic- 
ular p± to the interface (j p ). We find the surprising re- 
sult, that "f w depends on the functional form chosen to fit 
the order parameter profile through the interface. In par- 
ticular, fitting the profile to an error function gives results 
for j w which are in excellent agreement with j p . How- 
ever, fitting our data to tanh(2z/u>), a functional form 
derived from mean- field arguments H, gives results for 
7u, which are systematically 15% smaller than j p . Since 
the tanh function is often used to fit interfacial profiles 
at the liquid/ vapor interface |13 19(| , this difference is im- 
portant to understand. 

For this study we perform continuous-space, molecular 
dynamics simulations on a system of particles interacting 
through a standard (12-6) Lennard- Jones potential. The 
potential between particles i and j takes the form 
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where ry is the distance between particles i and j, and 
e and a set the energy and length scales of the potential 
respectively. Here we take a cut-off of r c =2.5cx. Increas- 
ing r c merely shifts T c to higher values, which should 
have little effect on the capillary-wave properties while 
increasing computation time significantly. The trajecto- 
ries of the N particles of mass m, are obtained by step- 
wise integration of Newton's equations of motion (EOM) 
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T [e/ka] 


N 


L[a] 


L ± [a] 


time [t] 


0.8 


7200 


12.8 


127.0 


6000 




24000 


24.7 


127.0 


6000 




69360 


42.0 


127.0 


6000 




154400 


54.5 


195.6 


5000 




506880 


94.0 


125.6 


2800 


0.9 


14400 


15.1 


216.1 


11000 




40000 


25.2 


216.1 


10000 




115660 


42.9 


216.1 


7800 




154400 


54.5 


195.6 


5800 




506880 


94.0 


144.4 


4200 




1240000 


134.6 


164.4 


4200 


1.0 


14400 


13.3 


264.2 


16500 




48000 


25.7 


264.2 


14400 




138720 


43.7 


264.2 


10500 




170000 


54.5 


195.6 


13500 




590000 


94.0 


293.9 


4900 



TABLE I. Values used for the parameters of simulations: 
temperature T, number of particles N, L = Lu, Lj_, and 
duration of run. 




FIG. 1. Typical configuration of an equilibrated liq- 
uid/vapor interface at T=0.8 e/fca- Length of square cross 
section holding the interface is L— 12.8a. 



In addition to the force derived from the LJ potential, 
the EOM contains a velocity-dependent damping term 
and a noise term representing a viscous drag force and a 
weak stochastic force, respectively. The noise term W(t) 
is taken from a uniform distribution, which mean value 
is set from the temperature T and the damping coeffi- 
cient r through the fluctuation-dissipation theorem [5o| . 
The combination of the viscous damping and stochas- 
tic force terms in the EOM effectively couples the sys- 
tem to a heat bath. Our simulations are performed in 
the canonical ensemble with fixed particle number and 
volume (constant-NVT). The EOM for each particle is 
integrated with the velocity- Verlet [|T) algorithm with 
a time step At = 0.006r, where T—a(m/ 'e) 1 / 2 fixes the 
time scale. We set T=0.5r~ 1 . All results presented here 
are measured in reduced units, as derived from the fun- 
damental scales fixed by a, e, m, and the Boltzmann 
constant ks ■ To reduce computation time we use a com- 
bination of the Verlet and linked-cell list algorithms [pi] . 

Periodic boundary conditions are used in all three 
space dimensions, thus forcing the creation of (at least) 
two interfaces in a two-phase system. The system sizes, 
temperatures, particle numbers and equilibration times 
of our simulations are listed in Table [j . In Table | and 
throughout this paper, L refers to the dimensions of the 
square cross section parallel to the interface, which lies 
in the xy plane. Thus, L x =L y =L\\=L. Lj_ refers to 
the dimension of the box perpendicular to the plane of 
the interface. Simulations are performed for L ranging 
from 12.8cr to 134. Qa. The largest system we could run 



T (e/ka) 


Pv 


PL 


7u>e 


7p 


0.8 


0.020(1) 


0.730(1) 


0.37(3) 


0.39(1) 


0.9 


0.045(1) 


0.663(1) 


0.22(1) 


0.22(1) 


0.95 


0.066(1) 


0.623(1) 




0.15(1) 


1.0 


0.098(2) 


0.571(1) 


0.097(2) 


0.08(1) 



TABLE II. Calculated values of the bulk densities and 
surface tensions for different simulation temperatures. 



contains 1.24 million particles. After that, computation 
becomes prohibitively slow due to the large number of 
particles. At the other end of our size range, systems 
with L < 12<T demonstrate non-negligible finite-size ef- 
fects. Figure [l] shows a typical configuration of an equi- 
librated system of L=12.8cr at T=0.8e/£;B. 

Initial systems were built as follows: for each system 
size and temperature, we construct a slab of the liquid 
phase and center it in the middle of the simulation cell 
with the interface perpendicular to the z direction. L±_ 
is set such that it is at least twice the length of the 
liquid slab, allowing sufficient space for the bulk liquid 
and vapor densities to achieve constant values. Since the 
phase coexistence diagram is well known for this system 
[jlT| , P~8 22 1, we adjusted the density of the liquid and va- 
por regions to be close to their reported values for each 
temperature T. 
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FIG. 2. Averaged density profile after equilibration for 
T=0.8 and 1.0 e/k B for L=41.9a. 




FIG. 3. Coexistence curve: bulk density values are 
obtained from tail values of interfacial liquid/vapor den- 
sity profiles. Curve is a best fit using the following ex- 
pressions El, 0.5(p L + /9v)cr 3 =0.544 - 0.210k B T/e and 
(Pl - pv)° J =A(l - (T/r c )) - 318 . The best fit parameters 
are A=1.07 and T c =1.085. 



After the system has equilibrated the density profile 
is measured, i.e. the xy-cross-section averaged number 
density p{z) as a function of z. Over the course of a 
simulation for a given T and L, p[z) is measured every 
400 time steps. Once the interfaces have equilibrated, 
the density profiles are averaged over 10 5 - 2 x 10 5 Ai. 
Great care must be exercised in the averaging procedure. 
For each profile, the position of the interface is located 
to insure that the averaging does not artificially broaden 
the interface width due to drift in the interface positions. 
Figure || shows an example of an equilibrated, averaged 
density profile for T = 0.8 and l.Oe/fcs for L=41.9a. 

Bulk values for the density are extracted from the tail 
values (obtained through a fit described below) of the in- 
terfacial density profiles. Our final equilibrated values for 
the bulk liquid and vapor densities are listed in Table 0, 
and agree very well with values reported in the litera- 
ture. The derived coexistence curve, along with a fit to 
an expression suggested in Ref. |lj] are shown in Fig. [|. 
The very good agreement of the bulk values suggests that 
our systems are well equilibrated that our measurement 
procedures are sound. Since the simulations are started 
near their respective liquid and vapor values, the bulk 
density values shown in Fig. || all attained equilibrium 
values quickly. However, the interface structure did not 



equilibrate until the simulations had been run for the 
much longer times shown in Table ||. 

An important quantity characterizing the interface is 
the width. The intrinsic width of an interface is due 
to the intermixing of the two phases, which always oc- 
curs to a certain degree at finite, subcritical tempera- 
tures. In addition to this mixing, capillary-wave theory 
predicts that thermal fluctuations of the location of 
the interface will contribute to the total, cross-section av- 
eraged, measured width. This broadening depends pri- 
marily on the surface tension, the temperature, and the 
cross-sectional size of the interface, and the spatial di- 
mension. As an example, capillary-wave theory states 
that any two-dimensional crystal is unstable against ther- 
mal fluctuations (23) . 

Fluctuations in £(x,y), the mean location of the in- 
terface in the z direction, induces fluctuations in the to- 
tal area of the interface and can be easily determined 
by expanding the shape of the interface to first order. 
This approximation is accurate provided the interface is 
smooth, with no overhangs. The free energy of the inter- 
face is the product of its surface area and an interfacial 
energy density 7, which is assumed to be independent of 
local curvature. Fluctuations due to capillary waves have 
an energy cost due to the increase in the surface area of 
the interface. The resulting interfacial Hamiltonian can 
be expressed as the product of surface tension times the 
increase in interfacial area 



H{(} = 7 / dxdy 




j I dxdy \V((x,y)\ 2 
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(4) 



The capillary-wave spectrum can be calculated by sub- 
stituting the Fourier transform of C giving 



n{(} 



dqq 2 \((q)\ 2 



(5) 



where q represents a two-dimensional vector in reciprocal 
space, and T {((x, y)] =C(<1) is the Fourier transform of 
The equipartition theorem dictates the mean- 
square amplitude for each interfacial excitation mode, 
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and summing over all allowed modes, one gets 

k B T 
4tt 2 7 j% 



icr 



dq 

2 ' 

q 



2^7 \B 



(6) 

(7) 
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where q m \- a =2ir / L and <7 max =27r/i? Note that both lower 
and upper cut-offs are required to prevent the value of 
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the integral from diverging. The long-wavelength cutoff 
<7 m i n , is determined by Lu ^4j. The interpretation of the 
short- wavelength cutoff, q max , is not as clear. Werner et 
al. (ljjl^l have studied the dependence of B for poly- 
mer/polymer interfaces and suggest B scales inversely 
with the molecular weight. However, the exact nature of 
this short-wavelength cut-off remains an open question. 

In both simulations and experiments, the quantity 
measured is the total interfacial width, which includes 
contributions from the intrinsic width and the broad- 
ening due to capillary-wave fluctuations. The two ef- 
fects can be distinguished if one assumes that capillary- 
wave fluctuations are decoupled from the intrinsic pro- 
file. Therefore, the total interface profile ^f(z) may be 
expressed as a convolution of the intrinsic interface pro- 
file ip(z) and the effect due to capillary waves |fl6| , 



ip(z - z a ) P{z ) dz Q 



(9) 



Here, P(z Q ) is the probability of finding the interface at 



i.e. 



L X Ly JQ 



dx dy S(C(x,y) - z Q ). (10) 



The interfacial order parameter profile ^f(z) is related 
to the cross-section averaged density profile p(z) by the 
function 



9{z) 



PL - PV 



p{z) 



PL + PV 



(11) 



which scales the density profile so that *B(z) varies be- 
tween — 1 and 1. The variance of the derivative of the 
total profile d$>(z)/dz = <!>' can be used as a measure of 
the width of the interface. The variance of a distribution 
/ is given by 



v[f] 



JT2f(q)\q=0 



dq 



/(0) 



(12) 



where f(q) is the Fourier transform of f(z). Making use 
of the convolution theorem and Eqs. (^J) and ([l2|) it can 
be shown that |l(| 



vfif'] = v[ip'] 
A 2 = A 2 4 



+ v[P], 
— lnf — 



(13) 



The squared widths of the total and intrinsic interfacial 
profiles have been defined as A 2 = v[*f?'] and A 2 = v[ip'], 
respectively. Note that the average squared fluctuations 
of the interface about its mean location in the z direc- 
tion can be directly identified as (|£| 2 )=i;[P]. Thus, our 
choice of measure for the interfacial width clearly shows 
that the total interfacial width can be written as the sum 



of an intrinsic part and a contribution due to capillary- 
wave fluctuations. 

In order to verify this prediction, we performed sev- 
eral simulations on different system sizes. Traditionally, 
the order parameter interfacial profile has been fit with 
f(z) = tanh(2z/w t ) or an error function erf (y/nz/w e ). 
Using our data we can test these two fitting functions 
and the resulting predictions for 7. Another reason for 
fitting our results for ^f(z) to one of these two functions 
is that we found we can determine a value for A 2 more 
accurately than by extracting it directly from the data; 
once the fitting parameters of f(z) have been determined, 
v[f] can be easily calculated. The two different fitting 
functions we tested are 



f e (z,w e ) = erf 



2z 

ft(z,w t ) = tanh ( — 



(14) 



For these two functions, the variance of each in terms of 
their associated widths TOj's are 



A 



v[f' e ] = w 2 e /2n 
v[f' t ] = ^ 2 W 2 /48. 



(15) 



The simulations are performed for three temperatures, 
T=0.8, 0.9, and 1.0 e/fes. This range of temperatures 
is selected because at lower temperatures the interfacial 
width is comparable to the average interparticle distance, 
and therefore is difficult to measure accurately. The up- 
per bound is set by T c w 1.085e/fce. For each value of 
T, the profiles are fit to both f x 's described above. Near 
the interface, the fitting functions can hardly be distin- 
guished. In fact, some studies have used an error function 
for theoretical derivations while using a hyperbolic tan- 
gent function to fit their data 

We fit our interfacial profiles for data near the inter- 
face and data deep into the bulk liquid and vapor regions. 
There is no a priori requirement that a liquid/vapor den- 
sity profile must be symmetric about the center of the 
interface. However, we detected no significant amounts 
of asymmetry. For each temperature and system size, the 
simulations are run until the interfacial profiles show a 
constant A 2 . 

Figure |^ summarizes the analysis from our extensive 
molecular dynamics simulations of a liquid/ vapor inter- 
face. For both tanh and erf fits, the data confirm a loga- 
rithmic dependence of A 2 on system size. Using Eq. dl3| ) 
the surface tensions can be calculated from the slopes of 
the best fit lines in Fig. || ['y w =k-QT/ (2tt slope)] . The tem- 
perature dependence of the interfacial surface tensions 
calculated from our simulations are shown in Fig. ^. We 
compare these values of the surface tension with an in- 
dependent measurement obtained from the components 
of the pressure tensor, r y p =L±(p± —p\\) p5[, represented 
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FIG. 4. Surface tension vs temperature T. Results for the 
surface tension -y w obtained from fitting interfacial profiles 
with a tanh or erf are compared to the values calculated from 
the components of the pressure tensor, "/ p =L±(p± ~P\\). The 
error function results are in excellent agreement with 7 P while 
the tanh results are systematically 15% too low. The 7 P value 
for T=0.85 is taken from Ref. Q. The -y p value for T=0.95 
is calculated from a single simulation with L=41.9er. 
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In (L/cr) 

FIG. 5. Variance A 2 versus InL for T=0.8 (squares), 0.9 
(triangles) and 1.0 e/k-g. (circles). The open and solid sym- 
bols are obtained using hyperbolic tangent and error function 
fits to the interfacial profiles, respectively. Lines are linear 
least-squares fits to the error function data. 

by solid squares. The agreement between j p and "f w ob- 
tained from the error function fits is very good. Using the 
tanh fits we obtain surface tensions that are systemati- 
cally 15% lower than those from the error function fits, 
which follows from their larger slopes shown in Fig. |^. 
Thus we obtain the somewhat unsettling result that the 
value of the A 2 and hence j w depends on the form of the 
fitting function used to fit if>(z). 

To investigate the systematic discrepancy between the 
tanh and erf fits, we performed an unweighted fit a hy- 
perbolic tangent function to data generated with an error 
function. The results are shown in Fig. [| The integrand 
in the numerator of the variance is plotted vs z for both 
functions. The integral of each of these functions is pro- 
portional to A 2 . From Fig. |^ one can see that the tails in 
the integrand of the tanh function contribute more sig- 
nificantly than the error function, hence the larger mea- 
sured values of A 2 from the tanh fits. We conclude that 
the tails of interfacial profiles are better captured by fits 




-4 -2 2 4 

FIG. 6. Results for an unweighted fit of a tanh function 
to a error function for w e = l. 



to an error function. 

In this paper, we presented results of extensive molecu- 
lar dynamics simulations of liquid/vapor interfaces. Our 
data confirm the capillary-wave description of the inter- 
face structure between a Lennard- Jones liquid and its 
vapor phase. When measuring the interfacial width by 
using second moments of the interfacial profile deriva- 
tives, we can extract values for the surface tension that 
agree very well with calculations from the components 
of the pressure tensor. We have also shown that the 
more robust method of extracting the second moment is 
through fits to an error function, since using a hyperbolic 
tangent leads to systematic errors. 

The results presented here are for an isolated liq- 
uid/vapor interface. The width of the liquid and va- 
por regions were carefully chosen so that there was no 
interference between the two interfaces. An interesting 
extension of this work is to study the effect of a nearby 
substrate on an interface. The effect of a wall on the in- 
terface can be modeled by adding a potential energy term 
to the interface Hamiltonian, H[d], which depends on the 
distance d between the interface and the wall. This term 
is calculated by integrating the potential energy between 
the microscopic constituents of two macroscopic objects, 
i.e. the interface and a semi-infinite wall. For pure LJ in- 
teractions, the potential energy is proportional to A/d 2 , 
where A is the Hamaker constant p(J. The Hamaker 
constant contains information about the strength of the 
microscopic potential, geometrical factors, and macro- 
scopic properties of the wall. Since this additional term 
in the Hamiltonian is quadratic in d, TL can therefore 
be diagonalized by a Fourier transform and the deriva- 
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tion of the capillary-wave spectrum is similar to the one 
presented here. The effect of the substrate is to cut off 
the long wavelength capillary-wave fluctuations so that 
A 2 no longer depends on L\\ for small d. The interplay 
between L\\ and d is an interesting question for which 
computer simulations such as these can directly address. 

We thank Frank van Swol for helpful discussions. San- 
dia is a multiprogram laboratory operated by Sandia Cor- 
poration, a Lockheed Martin Company, for the United 
States Department of Energy under Contract DE-AC04- 
94AL85000. 



off, where Ap is the difference in bulk density between 
the two phases. 
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